##########################################################################################
### Replication Code: All Analysis                                             ###########
### Title: Interstate Conciliation in the Shadow of Deep Historical Grievances ###########
### Authors: Kai Quek and Jiaqian Ni                                           ###########
### Date: January 30, 2026                                                     ###########
##########################################################################################


# Set-up ------------
## Clean the R environment and set the working directory
rm(list = ls())
setwd("~/Desktop/Projects/Survey 2023/data analysis") # set your working directory here, which should also contain the survey data

## load the packages 
library(readr)
library(dplyr)
library(tidyverse)
library(gplots)
library(ggplot2)
library(Rmisc)


# China Experiment ------------

## experimental groups
# apology == 1: Japan apologizes + China accepts 
# apology == 2: Japan apologizes + China rejects 
# apology == 3: Japan doesn't apologize + China accepts 
# apology == 4: Japan doesn't apologize + China rejects
# apology == 5: Japan kneels (sharp apology) + China accepts
# apology == 6: Japan kneels (sharp apology) + China rejects 

## load the cleaned dataset 
apchina <- readRDS("Conciliate_China.RData")
apchina$apology <- factor(apchina$apology, levels=c(1:6))


## Variables (China) ----------
## treatments (China)
### direct vs quasi apology 
apchina$direct <- NA
apchina$direct[apchina$apology==3 | apchina$apology==4] <- 0 # quasi-apology 
apchina$direct[apchina$apology==1 | apchina$apology==2] <- 1 # direct-apology
apchina$direct[apchina$apology==5 | apchina$apology==6] <- 2 # sharp-apology 
apchina$Direct <- factor(apchina$direct, levels=c(0,1,2), 
                         labels=c("Quasi", "Direct", "Sharp"))
table(apchina$Direct)


### positive vs negative Chinese response 
apchina$china <- NA
apchina$china[apchina$apology==1|apchina$apology==3|apchina$apology==5] <- 1 
apchina$china[apchina$apology==2|apchina$apology==4|apchina$apology==6] <- 0 
apchina$China <- factor(apchina$china, levels=c(0,1), 
                        labels=c("Negative", "Positive"))
table(apchina$China)


## DV measures (China, stage 1)

### perceived sincerity 
summary(apchina$apo2) 
apchina$sincere <- apchina$apo2
apchina$Sincere <- NA
apchina$Sincere[apchina$sincere > 3] <- 1
apchina$Sincere[apchina$sincere <= 3] <- 0 

### willingness to accept (s1)
summary(apchina$apo3a) 
apchina$accept <- apchina$apo3a
apchina$Accept <- NA
apchina$Accept[apchina$accept > 4] <- 1
apchina$Accept[apchina$accept <= 4] <- 0

## DV measures (stage 2)
### willingness to accept (s2)
summary(apchina$apo3b) 
apchina$accept2 <- apchina$apo3b
apchina$Accept2 <-NA
apchina$Accept2[apchina$accept2 > 4] <- 1
apchina$Accept2[apchina$accept2 <= 4] <- 0

apchina$accept_diff <- apchina$accept2 - apchina$accept # change in acceptance

### approval 
summary(apchina$apo1b) # disapproval [1=strongly approve, 7=strongly disapprove]
apchina$approve <- NA # reverse coding to reflect approval 
apchina$approve[apchina$apo1b==7] <- 1
apchina$approve[apchina$apo1b==6] <- 2
apchina$approve[apchina$apo1b==5] <- 3
apchina$approve[apchina$apo1b==4] <- 4
apchina$approve[apchina$apo1b==3] <- 5
apchina$approve[apchina$apo1b==2] <- 6
apchina$approve[apchina$apo1b==1] <- 7

apchina$Approve <- NA
apchina$Approve[apchina$approve > 4] <- 1 
apchina$Approve[apchina$approve <=4] <- 0 

### peaceful coexistence 
summary(apchina$apo4_1)
apchina$peace <- apchina$apo4_1
apchina$Peace <- NA
apchina$Peace[apchina$peace >5] <- 1
apchina$Peace[apchina$peace <=5] <- 0

### feelings toward Japan
summary(apchina$apo5_12) 
apchina$feeling <- apchina$apo5_12

### warm feelings toward Japan (0-100 feeling thermometer)
apchina$Warm <- NA
apchina$Warm[apchina$feeling > 50] <- 1
apchina$Warm[apchina$feeling <=50] <- 0 

## Statistics (China) ----------
### Summary of Results (China survey)
res_sincere <- summarySE(data=apchina, measurevar="Sincere", 
                         groupvars="Direct", na.rm = TRUE) #percentage sincerity by apology type
res_accept <- summarySE(data=apchina, measurevar="Accept", 
                        groupvars="Direct", na.rm = TRUE) #percentage acceptance (stage-1) by apology type
res_accept2 <- summarySE(data=apchina, measurevar="accept2", 
                         groupvars="apology", na.rm = TRUE) #stage-2 acceptance by group
res_accept_diff <- summarySE(data=apchina, measurevar="accept_diff", 
                             groupvars="apology", na.rm = TRUE) #within-subject change in acceptance by group
res_approve <- summarySE(data=apchina, measurevar="approve", 
                         groupvars="apology", na.rm = TRUE)
res_peace <- summarySE(data=apchina, measurevar="Peace", 
                       groupvars="apology", na.rm = TRUE) #percentage believing in peaceful relations by group
res_feeling <- summarySE(data=apchina, measurevar="feeling", 
                         groupvars="apology", na.rm = TRUE) #raw feeling scores by group
res_warm <- summarySE(data=apchina, measurevar="Warm", 
                      groupvars="apology", na.rm = TRUE) #percentage having warm feelings by group


### Hypothesis testing (China) 

## comparing perceived sincerity by apology type 
prop.test(
  x = c(sum(apchina$Sincere[apchina$direct == 0]),
        sum(apchina$Sincere[apchina$direct == 1])),
  n = c(sum(apchina$direct == 0),
        sum(apchina$direct == 1)),
  correct = FALSE) # (Direct vs Quasi)
prop.test(
  x = c(sum(apchina$Sincere[apchina$direct == 1]),
        sum(apchina$Sincere[apchina$direct == 2])),
  n = c(sum(apchina$direct == 1),
        sum(apchina$direct == 2)),
  correct = FALSE)# (Direct vs Sharp) # footnote 15 


##compare willingness to accept (%) 
prop.test(
  x = c(sum(apchina$Accept[apchina$direct == 0]),
        sum(apchina$Accept[apchina$direct == 1])),
  n = c(sum(apchina$direct == 0),
        sum(apchina$direct == 1)),
  correct = FALSE)# (stage 1, Direct vs Quasi)

prop.test(
  x = c(sum(apchina$Accept[apchina$direct == 1]),
        sum(apchina$Accept[apchina$direct == 2])),
  n = c(sum(apchina$direct == 1),
        sum(apchina$direct == 2)),
  correct = FALSE)# (stage 1, Direct vs Sharp)


# footnote 20 (paired t-test)
t.test(apchina$accept[apchina$china==1], 
       apchina$accept2[apchina$china==1], 
       alternative = "two.sided", mu=0, paired=TRUE, 
       var.equal=FALSE, conf.level=0.95)

# footnote 21 (paired t-test)
t.test(apchina$accept[apchina$china==0], 
       apchina$accept2[apchina$china==0], 
       alternative = "two.sided", mu=0, paired=TRUE, 
       var.equal=FALSE, conf.level=0.95)


## Compare approval 
## diff in approval [direct] when China responds positively vs. negatively
prop.test(
  x = c(sum(apchina$Approve[apchina$apology == 1]),
        sum(apchina$Approve[apchina$apology == 2])),
  n = c(sum(apchina$apology == 1),
        sum(apchina$apology == 2)),
  correct = FALSE)

## diff in approval [quasi] when China responds positively vs. negatively
prop.test(
  x = c(sum(apchina$Approve[apchina$apology == 3]),
        sum(apchina$Approve[apchina$apology == 4])),
  n = c(sum(apchina$apology == 3),
        sum(apchina$apology ==4)),
  correct = FALSE)

## diff in approval [sharp] when China responds positively vs. negatively
prop.test(
  x = c(sum(apchina$Approve[apchina$apology == 5]),
        sum(apchina$Approve[apchina$apology == 6])),
  n = c(sum(apchina$apology == 5),
        sum(apchina$apology == 6)),
  correct = FALSE)

## diff in approval (condition 1 vs. 3)
prop.test(
  x = c(sum(apchina$Approve[apchina$apology == 1]),
        sum(apchina$Approve[apchina$apology == 3])),
  n = c(sum(apchina$apology == 1),
        sum(apchina$apology == 3)),
  correct = FALSE)

## diff in approval (condition 2 vs. 4)
prop.test(
  x = c(sum(apchina$Approve[apchina$apology == 2]),
        sum(apchina$Approve[apchina$apology == 4])),
  n = c(sum(apchina$apology == 2),
        sum(apchina$apology == 4)),
  correct = FALSE)

# JAPAN Experiment ------------

## experimental groups 
# apology == 1: Japan apologizes + China accepts 
# apology == 2: Japan apologizes + China rejects 
# apology == 3: Japan doesn't apologize + China accepts 
# apology == 4: Japan doesn't apologize + China rejects
# apology == 5: Japan kneels (sharp apology) + China accepts
# apology == 6: Japan kneels (sharp apology) + China rejects 


## load the cleaned dataset 
apjapan <- readRDS("Conciliate_Japan.RData")
apjapan$apology <- factor(apjapan$apology, levels=c(1:6))

## Variables (Japan) ----------
## treatments (JAPAN)
### direct vs quasi apology
apjapan$direct <- NA
apjapan$direct[apjapan$apology==3 | apjapan$apology==4] <- 0 # quasi-apology
apjapan$direct[apjapan$apology==1 | apjapan$apology==2] <- 1 # direct-apology
apjapan$direct[apjapan$apology==5 | apjapan$apology==6] <- 2 # sharp-apology 
apjapan$Direct <- factor(apjapan$direct, levels=c(0,1,2), 
                         labels=c("Quasi", "Direct", "Sharp"))

### positive vs negative Chinese response 
apjapan$china <- NA
apjapan$china[apjapan$apology==1|apjapan$apology==3|apjapan$apology==5] <- 1 
apjapan$china[apjapan$apology==2|apjapan$apology==4|apjapan$apology==6] <- 0 
apjapan$China <- factor(apjapan$china, levels=c(0,1), 
                        labels=c("Negative", "Positive"))
table(apjapan$China)


## DV measures (Japan, stage 1)
### stage 1 support (1=strongly support)
summary(apjapan$apo1a) 
apjapan$approve <- NA # reverse coding to reflect approval 
apjapan$approve[apjapan$apo1a==7] <- 1 
apjapan$approve[apjapan$apo1a==6] <- 2 
apjapan$approve[apjapan$apo1a==5] <- 3 
apjapan$approve[apjapan$apo1a==4] <- 4 
apjapan$approve[apjapan$apo1a==3] <- 5 
apjapan$approve[apjapan$apo1a==2] <- 6 
apjapan$approve[apjapan$apo1a==1] <- 7 
table(apjapan$apo1a)
table(apjapan$approve)

### stage 1 percentage approval % 
apjapan$Approve <- NA
apjapan$Approve[apjapan$approve > 4] <- 1 
apjapan$Approve[apjapan$approve <=4] <- 0 

### perceived sincerity/seriousness (1=very sincere)
summary(apjapan$apo2) 
apjapan$sincere <- NA # reverse coding to reflect sincerity 
apjapan$sincere[apjapan$apo2==5] <- 1
apjapan$sincere[apjapan$apo2==4] <- 2
apjapan$sincere[apjapan$apo2==3] <- 3
apjapan$sincere[apjapan$apo2==2] <- 4
apjapan$sincere[apjapan$apo2==1] <- 5
table(apjapan$apo2)
table(apjapan$sincere)

### perceived sincerity (percentage) % 
apjapan$Sincere <- NA
apjapan$Sincere[apjapan$sincere > 3] <- 1
apjapan$Sincere[apjapan$sincere <= 3] <- 0

### to what extent do you think china is willing to accept
table(apjapan$apo3) 
# 1=completely, 3=neither, 5=not at all, 6=be, 7=a little, 9=not so much, 10=not
apjapan$accept <- NA
apjapan$accept[apjapan$apo3==5] <- 1
apjapan$accept[apjapan$apo3==10] <- 2
apjapan$accept[apjapan$apo3==9] <- 3
apjapan$accept[apjapan$apo3==3] <- 4
apjapan$accept[apjapan$apo3==7] <- 5
apjapan$accept[apjapan$apo3==6] <- 6
apjapan$accept[apjapan$apo3==1] <- 7
table(apjapan$apo3)
table(apjapan$accept)

### % perceived Chinese willingness to accept 
apjapan$Accept <- NA
apjapan$Accept[apjapan$accept > 4 ] <- 1
apjapan$Accept[apjapan$accept <= 4 ] <- 0

### stage 2 support (1=strongly support)
summary(apjapan$apo1b) 
apjapan$approve2 <- NA # reverse coding to reflect approval 
apjapan$approve2[apjapan$apo1b==7] <- 1 
apjapan$approve2[apjapan$apo1b==6] <- 2 
apjapan$approve2[apjapan$apo1b==5] <- 3 
apjapan$approve2[apjapan$apo1b==4] <- 4 
apjapan$approve2[apjapan$apo1b==3] <- 5 
apjapan$approve2[apjapan$apo1b==2] <- 6 
apjapan$approve2[apjapan$apo1b==1] <- 7 
table(apjapan$apo1b)
table(apjapan$approve2)

### stage 2 percentage approval %
apjapan$Approve2 <- NA
apjapan$Approve2[apjapan$approve2 > 4] <- 1
apjapan$Approve2[apjapan$approve2 <= 4] <- 0

### within-subject changes in approval 
apjapan$approve_diff <- apjapan$approve2 - apjapan$approve 

### belief in peaceful coexistence
summary(apjapan$apo4_1) 
apjapan$peace <- apjapan$apo4_1

## percentage of belief in peaceful coexistence % 
apjapan$Peace <- NA
apjapan$Peace[apjapan$peace > 5] <- 1
apjapan$Peace[apjapan$peace <= 5] <- 0

## feelings toward China
summary(apjapan$apo5_1) 
apjapan$feeling <- apjapan$apo5_1

## percentage with warm feelings % 
apjapan$Warm <- NA
apjapan$Warm[apjapan$feeling > 50 ] <- 1
apjapan$Warm[apjapan$feeling <=50] <- 0 


## Statistics (Japan) ----------
### Summary of Results (Japan survey)
jap_approve <- summarySE(data=apjapan, measurevar="approve", 
                         groupvars="apology", na.rm = TRUE)
jap_approve2 <- summarySE(data=apjapan, measurevar="approve2", 
                          groupvars="apology", na.rm = TRUE)
jap_approve_diff <- summarySE(data=apjapan, measurevar="approve_diff", 
                              groupvars="apology", na.rm = TRUE)
jap_sincere <- summarySE(data=apjapan, measurevar="Sincere", 
                         groupvars="Direct", na.rm = TRUE)#percentage of sincerity belief by apology type
jap_accept <- summarySE(data=apjapan, measurevar="Accept", 
                        groupvars="Direct", na.rm = TRUE)#percentage of perceived Chinese willingness to accept by apology type
jap_peace <- summarySE(data=apjapan, measurevar="Peace", 
                       groupvars="apology", na.rm = TRUE)
jap_feeling <- summarySE(data=apjapan, measurevar="feeling", 
                         groupvars="apology", na.rm = TRUE)
jap_warm <- summarySE(data=apjapan, measurevar="Warm", 
                      groupvars="apology", na.rm = TRUE)


### Hypothesis testing (Japan) 

## comparing perceived sincerity by apology type 
prop.test(
  x = c(sum(apjapan$Sincere[apjapan$direct == 0]),
        sum(apjapan$Sincere[apjapan$direct == 1])),
  n = c(sum(apjapan$direct == 0),
        sum(apjapan$direct == 1)),
  correct = FALSE) # (Direct vs Quasi)
prop.test(
  x = c(sum(apjapan$Sincere[apjapan$direct == 1]),
        sum(apjapan$Sincere[apjapan$direct == 2])),
  n = c(sum(apjapan$direct == 1),
        sum(apjapan$direct == 2)),
  correct = FALSE) # (Direct vs Quasi)# (Direct vs Sharp) # footnote 15 


## comparing perceived Chinese willingness to accept 
prop.test(
  x = c(sum(apjapan$Accept[apjapan$direct == 0], na.rm=TRUE),
        sum(apjapan$Accept[apjapan$direct == 1])),
  n = c(sum(apjapan$direct == 0),
        sum(apjapan$direct == 1)),
  correct = FALSE) # (Direct vs Quasi)
prop.test(
  x = c(sum(apjapan$Accept[apjapan$direct == 1]),
        sum(apjapan$Accept[apjapan$direct == 2])),
  n = c(sum(apjapan$direct == 1),
        sum(apjapan$direct == 2)),
  correct = FALSE) # (Direct vs Sharp)

## compare approval between stages (within-subject effects of official Chinese response) 
## McNemar tests (tests on paired binary data)
## condition 1 (direct + positive) p = 0.0002522
apjapn_c1 <- apjapan[apjapan$apology == 1, ]
tab_c1 <- table(apjapn_c1$Approve, apjapn_c1$Approve2)
mcnemar.test(tab_c1)

## condition 2 (direct + negative) p = 0.002089
apjapn_c2 <- apjapan[apjapan$apology == 2, ]
tab_c2 <- table(apjapn_c2$Approve, apjapn_c2$Approve2)
mcnemar.test(tab_c2)

## condition 3 (quasi + positive) p = 0.5419
apjapn_c3 <- apjapan[apjapan$apology == 3, ]
tab_c3 <- table(apjapn_c3$Approve, apjapn_c3$Approve2)
mcnemar.test(tab_c3)

## condition 4 (quasi + negative) p < 0.00001
apjapn_c4 <- apjapan[apjapan$apology == 4, ]
tab_c4 <- table(apjapn_c4$Approve, apjapn_c4$Approve2)
mcnemar.test(tab_c4)

## condition 5 (sharp + positive) p < 0.00001
apjapn_c5 <- apjapan[apjapan$apology == 5, ]
tab_c5 <- table(apjapn_c5$Approve, apjapn_c5$Approve2)
mcnemar.test(tab_c5)

## condition 6 (sharp + negative) p < 0.08199
apjapn_c6 <- apjapan[apjapan$apology == 6, ]
tab_c6 <- table(apjapn_c6$Approve, apjapn_c6$Approve2)
mcnemar.test(tab_c6)


# Sender-receiver gaps ------------
library(plyr)

### sender-receiver gap in sincerity perception 
## quasi-apology condition 
prop.test(
  x = c(sum(apchina$Sincere[apchina$direct == 0]),
        sum(apjapan$Sincere[apjapan$direct == 0])),
  n = c(sum(apchina$direct == 0),
        sum(apjapan$direct == 0)),
  correct = FALSE) 

## direct-apology condition 
prop.test(
  x = c(sum(apchina$Sincere[apchina$direct == 1]),
        sum(apjapan$Sincere[apjapan$direct == 1])),
  n = c(sum(apchina$direct == 1),
        sum(apjapan$direct == 1)),
  correct = FALSE) 

## sharp-apology condition 
prop.test(
  x = c(sum(apchina$Sincere[apchina$direct == 2]),
        sum(apjapan$Sincere[apjapan$direct == 2])),
  n = c(sum(apchina$direct == 2),
        sum(apjapan$direct == 2)),
  correct = FALSE) 


### sender-receiver gap in willingness to accept 
## quasi-apology
prop.test(
  x = c(sum(apchina$Accept[apchina$direct == 0]),
        sum(apjapan$Accept[apjapan$direct == 0], na.rm=TRUE)),
  n = c(sum(apchina$direct == 0),
        sum(apjapan$direct == 0)),
  correct = FALSE) 

## direct-apology
prop.test(
  x = c(sum(apchina$Accept[apchina$direct == 1]),
        sum(apjapan$Accept[apjapan$direct ==1])),
  n = c(sum(apchina$direct == 1),
        sum(apjapan$direct == 1)),
  correct = FALSE) 

## sharp
prop.test(
  x = c(sum(apchina$Accept[apchina$direct == 2]),
        sum(apjapan$Accept[apjapan$direct == 2])),
  n = c(sum(apchina$direct == 2),
        sum(apjapan$direct == 2)),
  correct = FALSE) 



# FIGURES ---------- 

## labels 
labs_full <- c("Direct +\nPositive", "Direct +\nNegative", 
               "Quasi +\nPositive", "Quasi + \nNegative", 
               "Sharp + \nPositive", "Sharp + \nNegative")
labs <- c("Quasi-apology", "Direct-apology", "Sharp-apology")


## FIGURE 3 ------
## Comparing Perceived Sincerity in China and Japan
china_sincere <- res_sincere %>% 
  mutate(upper = Sincere+ci, lower = Sincere-ci) %>% 
  dplyr::select(-c(N, sd, se, ci)) # china sincerity

japan_sincere <- jap_sincere %>% 
  mutate(upper = Sincere+ci, lower = Sincere-ci) %>% 
  dplyr::select(-c(N, sd, se, ci)) # japan sincerity

sincere <- rbind (china_sincere, japan_sincere)
sincere <- sincere %>% mutate(country=c(rep("China",3), rep("Japan", 3)))

ggplot(data=sincere, aes(x=Direct, y=Sincere*100, fill=country)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=lower*100, ymax=upper*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  theme_classic()+ scale_fill_grey(start=1, end=0.8)+
  labs(y="Japan's Apology is Sincere (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  geom_text(aes(label=paste(round(100*Sincere, 1), "%", sep="")), 
            vjust=6, size=3, position = position_dodge(0.9))+
  theme(legend.position=c(0.2, 0.9),
        legend.background = element_rect(color="black"),
        legend.title=element_blank(),
        legend.text = element_text(size=10),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 10),
        axis.title = element_text(colour = "black", size = 10)) 


## FIGURE 4 ------
## Actual and Perceived Chinese Willingness to Accept Apology
china_accept <- res_accept %>% 
  mutate(upper = Accept+ci, lower = Accept-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

japan_accept <- jap_accept %>% 
  mutate(upper = Accept+ci, lower = Accept-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

accept <- rbind (china_accept, japan_accept)
accept <- accept %>% mutate(country=c(rep("China",3), rep("Japan", 3)))

ggplot(data=accept, aes(x=Direct, y=Accept*100, fill=country)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=lower*100, ymax=upper*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  theme_classic()+ scale_fill_grey(start=1, end=0.8)+
  labs(y="(China is) Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  geom_text(aes(label=paste(round(100*Accept, 1), "%", sep="")), 
            vjust=6, size=3, position = position_dodge(0.9))+
  theme(legend.position=c(0.15, 0.9),
        legend.background = element_rect(color="black"),
        legend.title=element_blank(),
        legend.text = element_text(size=10),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 10),
        axis.title = element_text(colour = "black", size = 10)) 

## FIGURE 5 ------
## Changes in Willingness to Accept (China) 
p_accept_diff <- ggplot(data=res_accept_diff, aes(x=apology, y=accept_diff)) + 
  geom_point(shape=16, size=3)+
  geom_errorbar(aes(x=apology, ymin=accept_diff-ci, ymax=accept_diff+ci),
                width=0, position= position_dodge(width=0.8))+
  theme_classic()+
  scale_fill_grey(start=0.8, end=0.3)+
  labs(y="Changes in Willingness to Accept", x=NULL)+
  scale_x_discrete(labels=labs_full)+
  geom_hline(yintercept = 0, linetype = 2)
p_accept_diff


## FIGURE 6 ------
## Percentage Approval of China’s Response (China)
perc_approve <- summarySE(apchina, measurevar = "Approve", 
                          groupvars = "apology", na.rm = TRUE)
perc_approve <- perc_approve %>% mutate(china = rep(c("positive", "negative"), 3))
perc_approve$china <- as.factor(perc_approve$china)
p_perc_approve <- ggplot(perc_approve, aes(x=apology, y=Approve*100, fill=china)) +
  geom_bar(stat="identity", position = position_dodge(), color="black")+
  scale_fill_manual(values=c("grey80", "white"), name="China's response")+
  geom_errorbar(aes(ymin=(Approve-ci)*100, ymax=(Approve+ci)*100), 
                width=.1, position= position_dodge(0.9))+
  geom_text(aes(label=paste(round(100*Approve, 1), "%", sep=""), vjust=5))+
  ylab("Percentage Approval (%)") + 
  xlab(NULL) + theme_classic() +
  theme(legend.title=element_text(size=10),
        legend.text = element_text(size=9),
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))+
  scale_x_discrete(labels=labs)
p_perc_approve


## FIGURE 7 ------
## Percentage Approval of Japan’s Apology (Japan) 
perc_approve_jap <- summarySE(apjapan, measurevar = "Approve", 
                              groupvars = "apology", na.rm = TRUE) # Stage1
perc_approve2_jap <- summarySE(apjapan, measurevar = "Approve2", 
                               groupvars = "apology", na.rm = TRUE) # Stage2

approve_s1 <- perc_approve_jap %>% 
  mutate(upper = Approve+ci, lower = Approve-ci) %>% 
  dplyr::select(-c(N, sd, se, ci)) # stage1 95% confidence intervals

approve_s2 <- perc_approve2_jap %>% 
  mutate(upper = Approve2+ci, lower = Approve2-ci) %>% 
  dplyr::select(-c(N, sd, se, ci)) %>% 
  dplyr::rename(Approve = Approve2) # stage2 95% confidence intervals

approve_both <- rbind (approve_s1, approve_s2)
approve_both <- approve_both %>% 
  mutate(stage=c(rep("Stage 1",6), rep("Stage 2", 6))) #merge stage 1 &2 

ggplot(data=approve_both, aes(x=apology, y=Approve*100, fill=stage)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=apology, ymin=lower*100, ymax=upper*100),
                width=0.1, position= position_dodge(width=0.8), size=0.3)+
  theme_classic()+ scale_fill_grey(start=1, end=0.8)+
  labs(y="Percentage Approval (%)", x=NULL)+
  geom_text(aes(label=paste(round(100*Approve, 1), "%", sep="")), 
            vjust=6, size=2.5, position = position_dodge(0.9))+
  scale_x_discrete(labels=labs_full)+
  theme(legend.position=c(0.85, 0.9),
        legend.background = element_rect(color="black"),
        legend.title=element_blank(),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 8),
        axis.title = element_text(colour = "black", size = 8)) 



# APPENDIX B ------------

## Appendix B, Figure 1, likelihood of peaceful relations ------
china_peace <- res_peace %>% 
  mutate(upper = Peace+ci, lower = Peace-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

japan_peace <- jap_peace %>% 
  mutate(upper = Peace+ci, lower = Peace-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

peace <- rbind (china_peace, japan_peace)
peace <- peace %>% mutate(country=c(rep("China",6), rep("Japan", 6)))

## plot (china + japan) [Appendix B. Figure 1]
ggplot(data=peace, aes(x=apology, y=Peace*100, fill=country)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=apology, ymin=lower*100, ymax=upper*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  theme_classic()+ scale_fill_grey(start=1, end=0.8)+
  labs(y="China and Japan can Peacefully Co-exist and Co-prosper\n(% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs_full)+
  geom_text(aes(label=paste(round(100*Peace, 1), "%", sep="")), 
            vjust=6, size=3, position = position_dodge(0.9))+
  theme(legend.position=c(0.88, 0.9),
        legend.background = element_rect(color="black"),
        legend.title=element_blank(),
        legend.text = element_text(size=10),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 10),
        axis.title = element_text(colour = "black", size = 10)) 

## Appendix B, Figure 2, Feelings ------
china_feel <- res_warm %>% 
  mutate(upper = Warm+ci, lower = Warm-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

japan_feel <- jap_warm %>% 
  mutate(upper = Warm+ci, lower = Warm-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

feel <- rbind (china_feel, japan_feel)
feel <- feel %>% mutate(Country=c(rep("China",6), rep("Japan", 6)))

## plot (china + japan) [Appendix B. Figure 2]
ggplot(data=feel, aes(x=apology, y=Warm*100, fill=Country)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=apology, ymin=lower*100, ymax=upper*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  theme_classic()+ scale_fill_grey(start=1, end=0.8)+
  labs(y="Warm Feelings toward the Other Country (>50)\n(% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs_full)+
  geom_text(aes(label=paste(round(100*Warm, 1), "%", sep="")), 
            vjust=6, size=3, position = position_dodge(0.9))+
  theme(legend.position=c(0.88, 0.9),
        legend.background = element_rect(color="black"),
        legend.title=element_blank(),
        legend.text = element_text(size=10),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 10),
        axis.title = element_text(colour = "black", size = 10)) 



# SUBSET ANALYSIS ------------

## Subset analysis - China ------------
### Gender (1=male)
apchina <- apchina %>% mutate(Male = case_when(gender == 1 ~1, gender == 2 ~0))

### Education
apchina <- apchina %>% 
  mutate(University = case_when(edu1<=5 ~ 0, edu1>5 ~1))

### Political knowledge
# know 1: China's 2022 GDP growth (key=3%, know1==2)
apchina <- apchina %>% mutate(Poliknow1 = case_when(know1==2~1, TRUE~ 0))
# know 2: Governor of People's Bank of China (know2==1)
apchina <- apchina %>% mutate(Poliknow2 = case_when(know2==1~1, TRUE~ 0))
# know 3: Qingang's office, MFA minister (know3==3)
apchina <- apchina %>% mutate(Poliknow3 = case_when(know3==3~1, TRUE~ 0))
# know 4: Lee Kaun Yew (know4==3)
apchina <- apchina %>% mutate(Poliknow4 = case_when(know4==3~1, TRUE~ 0))

apchina <- apchina %>% 
  mutate(Poliknow = Poliknow1 + Poliknow2 + Poliknow3 + Poliknow4)


### nationalism 
# Everyone should support their country even when it is wrong
apchina <- apchina %>%mutate(nationalism=case_when(natwrong==7 ~ 1, 
                               TRUE ~ natwrong))
### hawkishness 
# Sometimes war is the only solution to international problems.
table(apchina$hawk)

### trust
table(apchina$trust)

### cosmopolitanism
# We are not just the citizens of our own countries, but also global citizens. 
table(apchina$cosmo_id)
apchina <- apchina %>% filter(cosmo_id>0)

### social dominance orientation
#(1) if everyone knows their position, many problems could be avoided.
apchina <- apchina %>% filter(sdo_1>0) %>% mutate(SDO_1=case_when(sdo_1>3~1, TRUE~0))
#(2) some people deserve better treatment than others
apchina <- apchina %>% filter(sdo_2>0) %>% mutate(SDO_2=case_when(sdo_2>3~1,TRUE~0))
#(3) some people are at the top, others bottom, this is necessary
apchina <- apchina %>% filter(sdo_3>0) %>% mutate(SDO_3=case_when(sdo_3>3~1, TRUE~0))
#(4) equality for everyone is good. 
apchina <- apchina %>% filter(sdo_4>0) %>% mutate(SDO_4=case_when(sdo_4<3~1, TRUE~0))

apchina <- apchina %>% mutate(sdo=SDO_1+SDO_2+SDO_3+SDO_4)

### China by age group ------
# 18-44 young; # 45 or above old
apchina <- apchina %>% mutate(young=case_when(age<45~1, TRUE~0))
accept_younger <- summarySE(data=apchina %>% filter(young==1), 
                            measurevar = "Accept2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         generation =rep("Younger", times=6))
accept_older<- summarySE(data=apchina %>% filter(young==0), 
                         measurevar = "Accept2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         generation =rep("Older", times=6))
accept_generation <- rbind(accept_younger, accept_older)

labs <- c("Direct", "Quasi", "Sharp")
ggplot(data=accept_generation, aes(x=Direct, y=Accept2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Accept2-ci)*100, ymax=(Accept2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(generation))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.20, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### China by nationalism ------
apchina <- apchina %>% 
  mutate(Nationalist = case_when(nationalism >3~1, TRUE~0))
accept_nationalist <- summarySE(data=apchina %>% filter(Nationalist==1), measurevar="Accept2", groupvars="apology", na.rm = TRUE)
accept_nationalist <- accept_nationalist %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         Nationalist = rep("Nationalist", times=6))
accept_non_nationalist <- summarySE(data=apchina %>% filter(Nationalist==0), measurevar="Accept2", groupvars="apology", na.rm = TRUE)
accept_non_nationalist <- accept_non_nationalist %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         Nationalist = rep("Non-Nationalist", times=6))
accept_nationalism <- rbind(accept_nationalist, accept_non_nationalist)

ggplot(data=accept_nationalism, aes(x=Direct, y=Accept2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Accept2-ci)*100, ymax=(Accept2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Nationalist))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  ylim(c(0,65))+
  theme(legend.position=c(0.80, 0.90),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### China by hawkishness ------
# "War is sometimes the only way to solve international problems." 
# 1 = strongly disagree, 5 = strongly agree
apchina <- apchina %>% mutate(hawkish=case_when(hawk>3~1, TRUE~0))

accept_hawkish <- summarySE(data=apchina %>% filter(hawkish==1), 
                            measurevar = "Accept2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Hawkish =rep("Hawks", times=6))
accept_dove <- summarySE(data=apchina %>% filter(hawkish==0), 
                         measurevar = "Accept2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Hawkish =rep("Doves", times=6))
accept_hawkdove <- rbind(accept_hawkish, accept_dove)

ggplot(data=accept_hawkdove, aes(x=Direct, y=Accept2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Accept2-ci)*100, ymax=(Accept2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Hawkish))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.20, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### China by education ------
accept_college <- summarySE(data=apchina %>% filter(University==1), 
                            measurevar = "Accept2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         edu =rep("University or above", times=6))
accept_no_college<- summarySE(data=apchina %>% filter(University==0), 
                              measurevar = "Accept2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         edu =rep("No University", times=6))
accept_education <- rbind(accept_college, accept_no_college)

ggplot(data=accept_education, aes(x=Direct, y=Accept2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Accept2-ci)*100, ymax=(Accept2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(edu))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.20, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))

### China by political knowledge ------
apchina <- apchina %>% mutate(Poliknow_high=case_when(
  Poliknow == 4 ~ "High", TRUE~"Low"))

accept_know_high <- summarySE(data=apchina %>% filter(Poliknow_high=="High"), 
                              measurevar="Accept2", groupvars="apology", na.rm = TRUE)
accept_know_high <- accept_know_high %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         Poliknow = rep("High Poliknow", times=6))
accept_know_low <- summarySE(data=apchina %>% filter(Poliknow_high=="Low"), measurevar="Accept2", groupvars="apology", na.rm = TRUE)
accept_know_low <- accept_know_low %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         Poliknow = rep("Low Poliknow", times=6))
accept_know <- rbind(accept_know_high, accept_know_low)

ggplot(data=accept_know, aes(x=Direct, y=Accept2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Accept2-ci)*100, ymax=(Accept2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Poliknow))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.20, 0.90),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### China by social dominance orientation ------
apchina <-apchina %>% mutate(SDO=case_when(sdo>1~1, TRUE~0))
accept_sdo_high <- summarySE(data=apchina %>% filter(SDO==1), 
                             measurevar="Accept2", groupvars="apology", na.rm = TRUE)
accept_sdo_high <- accept_sdo_high %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         SDO = rep("HIGH Social Dominance Orientation", times=6))
accept_sdo_low <- summarySE(data=apchina %>% filter(SDO==0), 
                            measurevar="Accept2", groupvars="apology", na.rm = TRUE)
accept_sdo_low <- accept_sdo_low %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         SDO = rep("LOW Social Dominance Orientation", times=6))
accept_sdo <- rbind(accept_sdo_high, accept_sdo_low)

ggplot(data=accept_sdo, aes(x=Direct, y=Accept2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Accept2-ci)*100, ymax=(Accept2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(SDO))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+ ylim(c(0,60))+
  theme(legend.position=c(0.80, 0.90),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### China by international trust ------
apchina <- apchina %>% mutate(trusting = case_when(trust >3 ~1, TRUE~0))
accept_trust <- summarySE(data=apchina %>% filter(trusting==1), 
                          measurevar="Accept2", groupvars="apology", na.rm = TRUE) %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         Trust = rep("High Trust", times=6))
accept_distrust <- summarySE(data=apchina %>% filter(trusting==0), 
                             measurevar="Accept2", groupvars="apology", na.rm = TRUE)
accept_distrust <- accept_distrust %>% 
  mutate(Direct = rep(c("Direct", "Indirect", "Sharp"), each=2), 
         China = rep(c("Positive", "Negative"), times=3), 
         Trust = rep("Low Trust", times=6))
accept_trusting <- rbind(accept_trust, accept_distrust)

ggplot(data=accept_trusting, aes(x=Direct, y=Accept2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Accept2-ci)*100, ymax=(Accept2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Trust))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Willing to Accept (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.80, 0.90),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))



## Subset analysis - Japan------------

### gender
apjapan <- apjapan %>% mutate(Male = case_when(gender == 1 ~1, gender == 2 ~0))

### education
table(apjapan$edu) # university = 5
apjapan <- apjapan %>% mutate(University = case_when(edu<5 ~ 0, edu>=5 ~1))

### nationalism 
# Everyone should support their country even when it is wrong
table(apjapan$natwrong)

### liberal/conservative (1= very liberal, 7= very conservative)
table(apjapan$leftright)

### hawkishness 
# Sometimes war is the only solution to international problems.
table(apjapan$hawk)

### cosmopolitan 
# We are not only citizens of Japan, we are citizens of the world.
# 1= strongly agree
table(apjapan$cosmo_id)
apjapan <- apjapan %>% 
  mutate(cosmo = case_when(
    cosmo_id==1~5, cosmo_id==2~4, cosmo_id==3~3, cosmo_id==4~2, cosmo_id==5~1))

### china is a big threat to japan
# 1= strongly disagree, 5= strongly agree
table(apjapan$ch_threat1)


### Japan by age group ------
# 18-44 young; 45 or above = old
apjapan <- apjapan %>% mutate(young=case_when(age<45~1, TRUE~0))
approve_younger <- summarySE(data=apjapan %>% filter(young==1), 
                             measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         generation =rep("Younger", times=6))
approve_older<- summarySE(data=apjapan %>% filter(young==0), 
                          measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         generation =rep("Older", times=6))
approve_generation <- rbind(approve_younger, approve_older)

ggplot(data=approve_generation, aes(x=Direct, y=Approve2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Approve2-ci)*100, ymax=(Approve2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(generation))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Approval Rate (% points)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.50, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### Japan by nationalism ------
apjapan <- apjapan %>% mutate(nationalistic=case_when(natwrong>2~1, TRUE~0))
approve_nationalist <- summarySE(data=apjapan %>% filter(nationalistic==1), 
                                 measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Nationalist =rep("Nationalist", times=6))
approve_non_nationalist <- summarySE(data=apjapan %>% filter(nationalistic==0), 
                                     measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Nationalist =rep("Non-nationalist", times=6))
approve_nationalism <- rbind(approve_nationalist, approve_non_nationalist)
ggplot(data=approve_nationalism, aes(x=Direct, y=Approve2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Approve2-ci)*100, ymax=(Approve2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Nationalist))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Approval Rate (% points)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.45, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### Japan by hawkishness ------
apjapan <- apjapan %>% mutate(hawkish=case_when(hawk>2~1, TRUE~0))
approve_hawkish <- summarySE(data=apjapan %>% filter(hawkish==1), 
                             measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Hawkish =rep("Hawks", times=6))
approve_dove <- summarySE(data=apjapan %>% filter(hawkish==0), 
                          measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Hawkish =rep("Doves", times=6))
approve_hawkdove <- rbind(approve_hawkish, approve_dove)
ggplot(data=approve_hawkdove, aes(x=Direct, y=Approve2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Approve2-ci)*100, ymax=(Approve2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Hawkish))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Approval Rate (% points)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.50, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### Japan by education ------
approve_college <- summarySE(data=apjapan %>% filter(University==1), 
                             measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         College =rep("University or above", times=6))
approve_non_college <- summarySE(data=apjapan %>% filter(University==0), 
                                 measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         College =rep("No University", times=6))
approve_edu <- rbind(approve_college, approve_non_college)
ggplot(data=approve_edu, aes(x=Direct, y=Approve2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Approve2-ci)*100, ymax=(Approve2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(College))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Approval Rate (% points)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.45, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))

### Japan by threat perception ------
apjapan <- apjapan %>% mutate(threat=case_when(ch_threat1>3~1, TRUE~0))
approve_threat_high <- summarySE(data=apjapan %>% filter(threat==1), 
                                 measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Threat =rep("High Threat Perception", times=6))
approve_threat_low <- summarySE(data=apjapan %>% filter(threat==0), 
                                measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Threat =rep("Low Threat Perception", times=6))
approve_threat <- rbind(approve_threat_high, approve_threat_low)
ggplot(data=approve_threat, aes(x=Direct, y=Approve2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Approve2-ci)*100, ymax=(Approve2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Threat))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Approval Rate (% points)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.20, 0.90),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### Japan by cosmopolitanism ------
apjapan <- apjapan %>% mutate(cosmopolitan=case_when(cosmo>3~1, TRUE~0))
approve_cosmo_high <- summarySE(data=apjapan %>% filter(cosmopolitan==1), 
                                measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Cosmo =rep("Cosmopolitan", times=6))
approve_cosmo_low <- summarySE(data=apjapan %>% filter(cosmopolitan==0), 
                               measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         Cosmo =rep("Non-Cosmopolitan", times=6))
approve_cosmo <- rbind(approve_cosmo_high, approve_cosmo_low)
ggplot(data=approve_cosmo, aes(x=Direct, y=Approve2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Approve2-ci)*100, ymax=(Approve2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(Cosmo))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Approval Rate (% points)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.80, 0.90),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


### Japan by conservatism ------
# 1= very liberal, 7 = very conservative
apjapan <- apjapan %>% mutate(conservative = case_when(leftright > 4 ~1, 
                                                       leftright < 4~0))
approve_conservative <- summarySE(data=apjapan %>% filter(conservative==1), 
                                  measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         leftright =rep("Conservative", times=6))
approve_liberal <- summarySE(data=apjapan %>% filter(conservative==0), 
                             measurevar = "Approve2", group="apology", na.rm = TRUE) %>% 
  mutate(Direct =rep(c("Direct", "Indirect", "Sharp"), each=2),
         China =rep(c("Positive", "Negative"), times=3),
         leftright =rep("Liberal", times=6))
approve_leftright <- rbind(approve_conservative, approve_liberal)

ggplot(data=approve_leftright, aes(x=Direct, y=Approve2*100, fill=China)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=(Approve2-ci)*100, ymax=(Approve2+ci)*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  facet_wrap(vars(leftright))+
  scale_fill_manual(values=c("grey90", "grey50")) + theme_classic()+
  labs(y="Approval Rate (% points)", x=NULL)+
  scale_x_discrete(labels=labs)+
  theme(legend.position=c(0.85, 0.95),
        legend.title=element_text(size=8),
        legend.text = element_text(size=8),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 9),
        axis.title = element_text(colour = "black", size = 9), 
        strip.background=element_rect(colour="black",
                                      fill=c("grey99")))


# USA Experiment ------------
## load the cleaned dataset
apusa <- readRDS("Conciliate_USA.RData")
apusa$apology <- factor(apusa$apology, levels=c(1:6))

## apology type
apusa$direct <- NA
apusa$direct[apusa$apology==3 | apusa$apology==4] <- 0 # quasi-apology 
apusa$direct[apusa$apology==1 | apusa$apology==2] <- 1 # direct-apology
apusa$direct[apusa$apology==5 | apusa$apology==6] <- 2 # sharp-apology
apusa$Direct <- factor(apusa$direct, levels=c(0,1,2), 
                       labels=c("Quasi", "Direct", "Sharp"))
table(apusa$Direct)

## China's official response 
apusa$china <- NA
apusa$china[apusa$apology==1|apusa$apology==3|apusa$apology==5] <- 1 
apusa$china[apusa$apology==2|apusa$apology==4|apusa$apology==6] <- 0 
apusa$China <- factor(apusa$china, levels=c(0,1), 
                      labels=c("Negative", "Positive"))
table(apusa$China)

## DV - Perceived sincerity 
apusa$sincere <- apusa$apo2
apusa$Sincere <- NA
apusa$Sincere[apusa$sincere >3 ] <- 1 
apusa$Sincere[apusa$sincere <=3] <- 0 

## DV - feelings 
summary(apusa$apo5_1) # feelings toward China
apusa$feeling_china <- apusa$apo5_1

summary(apusa$apo5_4) # feelings toward Japan
apusa$feeling_japan <- apusa$apo5_4

## Summary of results 
usa_sincere <- summarySE(data=apusa, measurevar="Sincere", 
                         groupvars="Direct", na.rm = TRUE)
res_feel_china <- summarySE(data=apusa, measurevar="feeling_china", 
                            groupvars="apology", na.rm = TRUE)
res_feel_jap <- summarySE(data=apusa, measurevar="feeling_japan", 
                          groupvars="apology", na.rm = TRUE)

### Figure 8 in Appendix F ------
# Perceived Sincerity in China, Japan and the United States
china_sincere <- res_sincere %>% 
  mutate(upper = Sincere+ci, lower = Sincere-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

japan_sincere <- jap_sincere %>% 
  mutate(upper = Sincere+ci, lower = Sincere-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

usa_sincere <- usa_sincere |> 
  mutate(upper = Sincere+ci, lower = Sincere-ci) %>% 
  dplyr::select(-c(N, sd, se, ci))

sincere_all <- rbind(china_sincere, japan_sincere, usa_sincere)
sincere_all <- sincere_all |> mutate(Country=c(rep("China",3), rep("Japan", 3), 
                                               rep("USA", 3)))

ggplot(data=sincere_all, aes(x=Direct, y=Sincere*100, fill=Country)) + 
  geom_bar(stat="identity", position = position_dodge(0.9), color="black")+
  geom_errorbar(aes(x=Direct, ymin=lower*100, ymax=upper*100),
                width=0.1, position= position_dodge(width=0.9), size=0.3)+
  theme_classic()+ scale_fill_grey(start=1, end=0.6)+
  labs(y="Japan's Apology is Sincere (% of respondents)", x=NULL)+
  scale_x_discrete(labels=labs)+
  geom_text(aes(label=paste(round(100*Sincere, 1), "%", sep="")), 
            vjust=6, size=3, position = position_dodge(0.9))+
  theme(legend.position=c(0.2, 0.95),
        legend.background = element_rect(color="black"),
        legend.title=element_blank(),
        legend.text = element_text(size=10),
        legend.direction = "horizontal", 
        legend.key.size = unit (0.4, 'cm'), 
        axis.text = element_text(colour = "black", size = 10),
        axis.title = element_text(colour = "black", size = 10)) 


### Figure 9 in Appendix F ------
# US Public’s Feelings toward China and Japan
res_feel_PRC <- res_feel_china %>% 
mutate(china_upper = feeling_china+ci, 
       china_lower = feeling_china-ci) %>% 
  dplyr::select(-c(N, sd, se, ci)) %>% 
  dplyr::rename(feel = feeling_china, upper=china_upper, lower=china_lower)

res_feel_JAP <- res_feel_jap %>% 
  mutate(jap_upper=feeling_japan+ci, 
         jap_lower=feeling_japan-ci) %>% 
  dplyr::select(-c(N, sd, se, ci)) %>% 
  dplyr::rename(feel=feeling_japan, upper=jap_upper, lower=jap_lower)

res_feel <- rbind (res_feel_PRC, res_feel_JAP)
res_feel <- res_feel %>% mutate(country=c(rep("China",6), rep("Japan", 6)))

ggplot(data=res_feel, aes(x=apology, y=feel, color=country)) + 
  geom_point(aes(shape=country), size=3, position=position_dodge(width=0.8))+
  geom_errorbar(aes(x=apology, ymin=lower, ymax=upper),
                width=0, position= position_dodge(width=0.8))+
  theme_classic()+
  labs(y="Feelings (0=cold, 100=warm)", x=NULL)+ylim(35,70)+
  scale_x_discrete(labels=labs_full)+
  geom_hline(yintercept = 50, linetype = 2)

### Hypothesis testing ------
# t tests on perceived sincerity 
prop.test(
  x = c(sum(apusa$Sincere[apusa$direct == 0], na.rm=TRUE),
        sum(apusa$Sincere[apusa$direct == 1], na.rm=TRUE)),
  n = c(sum(apusa$direct == 0),
        sum(apusa$direct == 1)),
  correct = FALSE) 
prop.test(
  x = c(sum(apusa$Sincere[apusa$direct == 1], na.rm=TRUE),
        sum(apusa$Sincere[apusa$direct == 2], na.rm=TRUE)),
  n = c(sum(apusa$direct == 1),
        sum(apusa$direct == 2)),
  correct = FALSE) 


# compare feelings toward Japan when China responds positively vs negatively 
t.test(apusa$feeling_japan[apusa$apology==1],
       apusa$feeling_japan[apusa$apology==2],
       alternative = "two.sided", mu=0, paired=FALSE,
       var.equal=FALSE, conf.level=0.95) # direct apology
t.test(apusa$feeling_japan[apusa$apology==3],
       apusa$feeling_japan[apusa$apology==4],
       alternative = "two.sided", mu=0, paired=FALSE,
       var.equal=FALSE, conf.level=0.95) # quasi apology
t.test(apusa$feeling_japan[apusa$apology==5],
       apusa$feeling_japan[apusa$apology==6],
       alternative = "two.sided", mu=0, paired=FALSE,
       var.equal=FALSE, conf.level=0.95) # sharp apology

# t test on feelings toward China
t.test(apusa$feeling_china[apusa$apology==1], 
       apusa$feeling_china[apusa$apology==2], 
       alternative = "two.sided", mu=0, paired=FALSE,
       var.equal=FALSE, conf.level=0.95) # direct apology (China positive vs negative)
t.test(apusa$feeling_china[apusa$apology==3], 
       apusa$feeling_china[apusa$apology==4], 
       alternative = "two.sided", mu=0, paired=FALSE,
       var.equal=FALSE, conf.level=0.95) # quasi apology (China positive vs negative)
t.test(apusa$feeling_china[apusa$apology==5], 
       apusa$feeling_china[apusa$apology==6], 
       alternative = "two.sided", mu=0, paired=FALSE,
       var.equal=FALSE, conf.level=0.95) # sharp apology (China positive vs negative)




# Structural Topic Models ------------
## Chinese open-ended responses -----------
rm(list = ls())
setwd("~/Desktop/Projects/Survey 2023/data analysis") # set your working directory here, which should also contain the open-ended response data

## Load the required packages
library(stm)       # version 1.3.6
library(tm)        # version 0.7-8
library(stringi)   # version 1.7.6
library(extrafont) # version 0.17
library(dplyr)

## Import the dataset that contains all valid open-ended responses
data <- read.csv("china_open.csv", na.strings = c("", "NA"))
data <- data %>% mutate(quality= case_when(direct==0~0, TRUE~1))

## Process the text data
processed <- textProcessor(data$translation, metadata = data)
out <- prepDocuments(processed$documents, processed$vocab, processed$meta, lower.thresh = 3)
docs <- out$documents
vocab <- out$vocab
meta <- out$meta

## Run the STM (estimating a five-topic STM)
set.seed(22)
selectedmodel <- stm(out$documents, out$vocab, K = 5,
                     prevalence =~ direct * china,
                     data = out$meta, init.type = "Spectral", 
                     verbose = FALSE)
plot(selectedmodel)
## Display the words that characterize each topic 
labelTopics(selectedmodel)

thoughts1 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 1)$docs[[1]] # magnanimity
thoughts2 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 2)$docs[[1]] # support country
thoughts3 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 3)$docs[[1]] # mixed
thoughts4 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 4)$docs[[1]] # history
thoughts5 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 5)$docs[[1]] # peace

prep <- estimateEffect(1:5 ~ quality * china + approve, selectedmodel,
                       meta = out$meta,
                       uncertainty = "None")

plot(prep, covariate="china", topics=c(1,2,3,4,5),
     model=selectedmodel, method="difference", 
     cov.value1 = "1", cov.value2 = "0", 
     xlab="Change in topical prevalence from negative to positive",
     main= "Effect of China's Reaction", labeltype="custom", 
     custom.labels = c("Chinese magnanimity", "support country", "mixed", 
                       "history", "peace & development" ))

plot(prep, covariate="quality", topics=c(1,2,3,4,5),
     model=selectedmodel, method="difference", 
     cov.value1 = "1", cov.value2 = "0", 
     xlab="Change in topical prevalence from quasi to direct/sharp",
     main= "Effect of Apology Quality", labeltype="custom", 
     custom.labels = c("Chinese magnanimity", "support country", "mixed", 
                       "history", "peace & development" ))

## Japanese open-ended responses ------------
rm(list = ls())
setwd("~/Desktop/Projects/Survey 2023/data analysis") # set your working directory here, which should also contain the open-ended response data

## load the dataset containing English translations
japan_open_translated <- read.csv("jap_open.csv", na.strings = c("", "NA")) %>% 
  mutate(quality= case_when(direct==0~0, TRUE~1))

## Process the text data
processed <- textProcessor(japan_open_translated$translation, 
                           metadata = japan_open_translated)
out <- prepDocuments(processed$documents, processed$vocab, 
                     processed$meta, lower.thresh = 2)
docs <- out$documents
vocab <- out$vocab
meta <- out$meta


## Run the STM (estimating a 6-topic model)
set.seed(42)
selectedmodel <- stm(out$documents, out$vocab, K = 6,
                     prevalence =~ quality * china + Approve2,
                     data = out$meta, init.type = "Spectral", 
                     verbose = FALSE)

plot(selectedmodel)

## Display the words that characterize each topic 
labelTopics(selectedmodel)

set.seed(42)
thoughts1 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 100, 
                          topics = 1)$docs[[1]]
thoughts2 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 2)$docs[[1]]
thoughts3 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 3)$docs[[1]]
thoughts4 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 4)$docs[[1]]
thoughts5 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 5)$docs[[1]]
thoughts6 <- findThoughts(selectedmodel, 
                          texts = meta$translation, 
                          n = 50, 
                          topics = 6)$docs[[1]]

## plot the results 
set.seed(42)
prep <- estimateEffect(1:6 ~ quality * china + Approve2, selectedmodel,
                       metadata = out$meta, uncertainty = "None")

plot(prep, covariate="china", topics=c(1,2,3,4,5,6),
     model=selectedmodel, method="difference", 
     cov.value1 = "1", cov.value2 = "0", 
     xlab="Change in topical prevalence from negative to positive",
     main= "Effect of China's Reaction", labeltype = "custom", width=10,
     custom.labels = c("Improve relations", "Apology fatigue", "No need to apologize", 
                       "Mixed", "Can't trust China", "Apology worked"))

plot(prep, covariate="quality", topics=c(1,2,3,4,5,6),
     model=selectedmodel, method="difference", 
     cov.value1 = "1", cov.value2 = "0", 
     xlab="Change in topical prevalence from quasi to direct/sharp",
     main= "Effect of Apology Quality", 
     width=10,labeltype = "custom", 
     custom.labels = c("Improve relations", "Apology fatigue", "No need to apologize", 
                       "Mixed", "Can't trust China", "Apology worked"))



